Alejandro Naranjo Caraza
Orlando Almazán

import random
import plotly.graph_objects as go
import cvxpy as cp

import matplotlib.pyplot as plt
import numpy as np
from scipy import linalg, optimize

from AW_loadProblem import loadProblem
from mActiveSet import *

Imported: loadProblem (by AW).

---

I. Problema general de programación cuadrática

Un problema general de programación cuadrática (QP, por sus siglas en inglés) en su forma estándar de minimización puede formularse como:

$$ \begin{aligned} \min_{x \in \mathbb{R}^n} \quad & \frac{1}{2} x^\top G x + c^\top x \\ \text{sujeto a} \quad & A_i x = b_i, \quad \forall i \in \varepsilon, \\ & A_i x \leq b_i, \quad \forall i \in \mathcal{I}, \end{aligned} $$

donde:

• $x \in \mathbb{R}^n$ es el vector de variables de decisión.
• $G \in \mathbb{R}^{n \times n}$ es una matriz simétrica definida positiva (o semidefinida positiva), que define el término cuadrático de la función objetivo.
• $c \in \mathbb{R}^n$ es el vector de coeficientes lineales en la función objetivo.
• $A \in \mathbb{R}^{m \times n}$ es la matriz de coeficientes de las restricciones.
• $b \in \mathbb{R}^m$ es el vector del lado derecho de las restricciones.
• $\varepsilon \subseteq \{1, \dots, m_1\}$ es el conjunto de índices correspondientes a restricciones de igualdad.
• $\mathcal{I} \subseteq \{1, \dots, m_2\}$ es el conjunto de índices correspondientes a restricciones de desigualdad.

Se asume que $\varepsilon \cap \mathcal{I} = \emptyset$ y que $\varepsilon \cup \mathcal{I} = \{1, \dots, m\}$. Es decir, cada restricción está clasificada de manera única como igualdad o desigualdad.

---

II. Implementación del problema

Definimos la clase QuadraticProgramming para representar un problema de programación cuadrática y evaluarlo, verificar su convexidad estricta y generar un punto factible.

__init__(...)

Inicializa el problema con:

• G_matrix: matriz cuadrática $G$.
• c_vec: vector lineal $c$.
• A_matrix, b_vec: matriz y vector de restricciones (opcionales).
• eps_dic, l_dic: índices de las restricciones de igualdad ($\varepsilon$) y desigualdad ($\mathcal{I}$).

evaluate(x)

Calcula el valor de la función objetivo en un punto $x$ dado:

$$ f(x) = \frac{1}{2} x^\top G x + c^\top x $$

generate_feasible()

Busca un punto factible resolviendo un problema lineal auxiliar con scipy.optimize.linprog.

• Si hay restricciones y están clasificados los índices (eps y l), separa las filas correspondientes.
• Si no hay restricciones, devuelve el vector cero.
• Si no encuentra solución, imprime un mensaje.

is_strict_convex()

Devuelve True si todos los autovalores de $G$ son positivos (es decir, $G$ es definida positiva), lo que implica que la función objetivo es estrictamente convexa:

class QuadraticProgramming:
    def __init__(
        self, G_matrix, c_vec, A_matrix=None, b_vec=None, eps_dic=None, l_dic=None
    ):
        self.G = G_matrix
        self.A = A_matrix
        self.c = c_vec
        self.b = b_vec
        self.restrictions_eps = eps_dic
        self.restrictions_l = l_dic

    def evaluate(self, x):
        res = (1 / 2) * x @ (self.G @ x) + self.c @ x
        return res

    def generate_feasible(self):
        if self.A is not None and self.b is not None:
            if self.restrictions_l or self.restrictions_eps:
                A_eps = self.A[self.restrictions_eps, :]
                A_l = self.A[self.restrictions_l, :]
                b_eps = self.b[self.restrictions_eps]
                b_l = self.b[self.restrictions_l]
                c = np.zeros(A_eps.shape[1])
                optim = optimize.linprog(c, A_ub=A_l, A_eq=A_eps, b_ub=b_l, b_eq=b_eps)
            else:
                optim = optimize.linprog(self.c, A_eq=self.A, b_eq=self.b)
            if optim.success:
                res = optim.x
            else:
                print("No feasible point found.")
                print(optim.message)
                res = None
        else:
            res = np.zeros(self.G.shape[1])
        return res

    def is_strict_convex(self):
        return bool(np.all(np.linalg.eigvals(self.G) > 0))

Función schur_complement(problem)

Esta función resuelve un problema cuadrático con restricciones de igualdad usando el método del complemento de Schur.

Pasos clave:

• Extrae la submatriz de restricciones de igualdad $A$ y el vector $b$ desde el objeto problem.
• Calcula la descomposición de Cholesky de $G$:
  $$ G = L^\top L $$
  
• Calcula $G^{-1}$ eficientemente usando $L$: $$ G^{-1} = (L^\top)^{-1} L^{-1} $$
• Define el sistema reducido: $$ M = L^{-1} A^\top $$
• Resuelve primero por los multiplicadores de Lagrange $\lambda$: $$ \lambda = -(M^\top M)^{-1} (b + A G^{-1} c) $$
• Luego obtiene $x$ óptimo: $$ x = -G^{-1} (c + A^\top \lambda) $$

Devuelve:

Una tupla (x, λ) con el punto óptimo y los multiplicadores.

Implementación

def schur_complement(problem):
    if problem.A is None or problem.b is None:
        raise ValueError("Schur complement expects constraints.")
    if problem.restrictions_eps or problem.restrictions_l:
        A = problem.A[problem.restrictions_eps, :]
        b = problem.b[problem.restrictions_eps]
    else:
        A = problem.A
        b = problem.b
    A_transp = np.transpose(A)
    G = problem.G
    c = problem.c
    L = np.transpose(linalg.cholesky(G))
    L_transp = np.transpose(L)
    L_inv = linalg.inv(L)
    M = linalg.solve(L, A_transp)
    M_transp = np.transpose(M)
    G_inv = linalg.inv(L_transp) @ L_inv
    lambd = linalg.solve(M_transp @ M, -b - A @ G_inv @ c)
    x = -G_inv @ (c + A_transp @ lambd)
    res = x, lambd
    return res

Función null_space(problem)

Resuelve el problema usando el método del espacio nulo, útil cuando las restricciones de igualdad son de rango completo.

Pasos clave:

• Verifica que $A$ tenga más columnas que filas y que tenga rango completo.

• Aplica descomposición QR a $A^\top$: $$ A^\top = Q R $$ donde:

  • $Y$ es la base del espacio columna de $A^\top$
  • $Z$ es la base del espacio nulo de $A$

• Descompone la solución en dos partes:

  • $x_p$: solución particular que satisface $Ax = b$
  • $x_h$: solución homogénea que vive en el espacio nulo de $A$
  $$ x = x_p + x_h $$

• Calcula: $$ x_p = Y R^{-\top} b \\ x_h = Z (Z^\top G Z)^{-1} Z^\top (-c - G x_p) \\ x = x_p + x_h $$

• Calcula los multiplicadores: $$ \lambda = R^{-1} (-Y^\top (c + G x)) $$

Devuelve:

Una tupla (x, λ) con la solución óptima y los multiplicadores. Not: Ambas funciones suponen que el problema tiene solo restricciones de igualdad y que $G$ es simétrica definida positiva.

Implementación

def null_space(problem):
    if problem.A is None or problem.b is None:
        raise ValueError("Null space expects constraints.")
    if problem.restrictions_eps or problem.restrictions_l:
        A = problem.A[problem.restrictions_eps, :]
        b = problem.b[problem.restrictions_eps]
    else:
        A = problem.A
        b = problem.b
    if A.shape[0] > A.shape[1]:
        raise ValueError("Null space error. Matrix shape.")
    rank_A = np.linalg.matrix_rank(A)
    if rank_A < A.shape[0]:
        raise ValueError("Null space error. Matrix not full row rank")
    A_transp = np.transpose(A)
    G = problem.G
    c = problem.c
    Q, R = linalg.qr(A_transp, mode="full")
    Y = Q[:, :rank_A]
    Z = Q[:, rank_A:]
    R = R[:rank_A, :]
    Z_transp = np.transpose(Z)
    R_transp = np.transpose(R)
    Y_transp = np.transpose(Y)
    cp = linalg.solve(R_transp, b)
    xp = Y @ cp
    ch = linalg.solve(Z_transp @ G @ Z, -Z_transp @ (c + G @ xp))
    x = xp + Z @ ch
    lambd = linalg.solve(R, -Y_transp @ (c + G @ x))
    res = x, lambd
    return res

Método del Conjunto Activo

Entradas:

• problem: instancia de QuadraticProgramming
• x0: punto inicial
• W0: conjunto activo inicial (índices de restricciones activas)
• tol: tolerancia de convergencia
• max_iter: número máximo de iteraciones
• trajectory: (opcional) si es True, guarda trayectoria de $x_k$ y $q_k$

Salidas:

• Solución óptima x*
• Conjunto activo final W*
• (Opcional) trayectoria de iteraciones

Estructura:

• active_set_branch_1: busca dirección $d_k$ resolviendo un subproblema QP reducido:

  • Usa Schur complement o espacio nulo según el caso.
  • Calcula paso $\alpha$ y actualiza $x_{k+1} = x_k + \alpha d_k$
  • Si alguna restricción desigual se activa, se añade a $W_k$.

• active_set_branch_2: revisa los multiplicadores de Lagrange $\lambda_k$:

  • Si alguno $\lambda_j < 0$ y $j \notin \varepsilon$, elimina $j$ de $W_k$
  • Si todos $\lambda_j \geq 0$, se ha encontrado solución.

Cicla entre rama 1 y rama 2 hasta converger o alcanzar max_iter.

Implementación

def active_set_branch_1(
    problem,
    x0,
    W0,
    zeros,
    x_trajectory=None,
    q_trajectory=None,
    tol=10e-10,
    max_iter=100,
):
    k = 0
    xk = x0
    Wk = W0
    next_iter = True
    if x_trajectory is not None:
        x_trajectory.append(xk)
    if q_trajectory is not None:
        q_trajectory.append(problem.evaluate(xk))
    while next_iter and k < max_iter:
        print(f"\n\033[4mBranch 1; k = {k}\033[0m")
        print(f"W{k} = {Wk}")
        gk = problem.G @ xk + problem.c
        if Wk:
            Ak = problem.A[Wk, :]
            zerosk = zeros[Wk]
            problem_k = QuadraticProgramming(problem.G, gk, Ak, zerosk)
            schur_condition_1 = problem_k.is_strict_convex()
            schur_condition_2 = (
                abs(problem_k.A.shape[1] - problem_k.A.shape[0])
                > problem_k.A.shape[1] / 2
            )
            if schur_condition_1 and schur_condition_2:
                dk, lambdk = schur_complement(problem_k)
                print("schur complement used")
            else:
                dk, lambdk = null_space(problem_k)
                print("null space used")
        else:
            problem_k = QuadraticProgramming(problem.G, gk)
            dk, lambdk = non_restriction_optimization(problem_k)
            print("non restriction optimization")
        dk_norm = np.linalg.norm(dk, ord=np.inf)
        print(f"norm(d{k}) = {np.round(dk_norm, decimals=3)}")
        if dk_norm > tol:
            Wk_set = set(Wk)
            temp_set_l = [i for i in problem.restrictions_l if i not in Wk_set]
            Al = problem.A[temp_set_l, :]
            bl = problem.b[temp_set_l]
            alpha_condition = (Al @ dk > 0) & (bl - Al @ xk > 0)
            temp_set_l = np.array(temp_set_l)[alpha_condition].tolist()
            if not temp_set_l:
                alpha = 1
                print("No alpha tilde.")

            else:
                Al = Al[alpha_condition, :]
                bl = np.extract(alpha_condition, bl)
                alpha_vec = (bl - Al @ xk) / (Al @ dk)
                index_min = np.argmin(alpha_vec)
                alpha_min = alpha_vec[index_min]
                j_min = int(temp_set_l[index_min])
                print(
                    f"alpha tilde = {
                        np.round(alpha_min, decimals=8)}"
                )
                alpha = min(alpha_min, 1.0)
            xk = xk + alpha * dk
            if x_trajectory is not None:
                x_trajectory.append(xk)
            if q_trajectory is not None:
                q_trajectory.append(problem.evaluate(xk))
            print(f"x{k+1} = {np.round(xk, decimals=3)}")
            print(f"q{k+1} = {np.round(problem.evaluate(xk), decimals=3)}")
            if alpha != 1:
                print(f"W{k+1} = W{k} U [{j_min}]")
                Wk.append(j_min)
                k += 1
            else:
                next_iter = False
                print(f"W{k+1} = W{k}")
        else:
            next_iter = False
    if next_iter and k == max_iter:
        raise MaxIterationsExceeded(
            f"Branch 1 maximum number of iterations ({max_iter}) exceeded."
        )
    else:
        res = xk, Wk, lambdk
    return res

def active_set_branch_2(problem, xk, Wk, lambdk):
    print("\n\033[4mBranch 2\033[0m")
    min_mu = 0
    min_j = None
    for i in range(len(Wk)):
        if Wk[i] not in problem.restrictions_eps and lambdk[i] < min_mu:
            min_mu = lambdk[i]
            min_j = Wk[i]
    if min_mu < 0:
        Wk.remove(min_j)
        next_iter = True
        print(f"m{min_j} = {min_mu} < 0. Back to Branch 1.")
    else:
        next_iter = False
        print("mj >= 0 forall j. Solution found.")
    return Wk, next_iter

def active_set(problem, x0, W0, tol=10e-9, max_iter=100, trajectory=False):
    if not isinstance(problem, QuadraticProgramming):
        raise TypeError("Expected Quadratic Programming Instance.")
    if not linalg.issymmetric(problem.G):
        raise ValueError("Optimization matrix is not symetric")
    if not set(problem.restrictions_eps).issubset(W0):
        raise ValueError("Error. W0 not subset of equality restrictions")
    if trajectory:
        x_trajectory = []
        q_trajectory = []
    else:
        x_trajectory = None
        q_trajectory = None
    xk = np.copy(x0)
    Wk = W0.copy()
    zeros = np.zeros(problem.A.shape[0])
    iter = 0
    next_iter = True
    while next_iter and iter < max_iter:
        print("\n-----------------------------------------------------------")
        print(f"ITERATION {iter}")
        print("-----------------------------------------------------------")
        xk, Wk, lambdk = active_set_branch_1(
            problem, xk, Wk, zeros, x_trajectory, q_trajectory
        )
        Wk, next_iter = active_set_branch_2(problem, xk, Wk, lambdk)
        iter += 1
    print("\n-----------------------------------------------------------")
    if next_iter and iter == max_iter:
        res = None
        print("Outer layer maximum iterations reached. No solution found.")
        print("-----------------------------------------------------------\n")
    else:
        if trajectory:
            x_trajectory = np.array(x_trajectory)
            q_trajectory = np.array(q_trajectory)
            res = xk, Wk, x_trajectory, q_trajectory
        else:
            res = xk, Wk
        print("ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND.")
        print("-----------------------------------------------------------\n")
        print(f"x* = {np.round(xk, decimals=3)}\nW* = {Wk}")
    return res

---

III. Problema 1

En el primer problema se pide minimizar la función objetivo:

$$ q(\vec{x}) = (x - 1)^2 + (y - 2.5)^2 $$

sujeta a las siguientes restricciones:

$$ \begin{aligned} -x + 2y - 2 &\leq 0 \\ x + 2y - 6 &\leq 0 \\ x - 2y - 2 &\leq 0 \\ -x &\leq 0 \\ -y &\leq 0 \end{aligned} $$

Primera prueba de optimización

Partimos del siguiente punto inicial:

$$ \vec{x}_0 = \begin{pmatrix} 2 \\ 0 \end{pmatrix} $$

y del conjunto activo inicial:

$$ W_0 = \{3\} \subset \{1, 2, 3, 4, 5\} $$

Nota: En la implementación del algoritmo, las restricciones están indexadas desde cero (0, 1, 2, ...), por lo que al ejecutar la función se toma $W_0 = \{2\}$.

G_mat = np.array([[2, 0], [0, 2]])
A_mat = np.array([[-1, 2], [1, 2], [1, -2], [-1, 0], [0, -1]])
c_vec = np.array([-2, -5])
b_vec = np.array([2, 6, 2, 0, 0])
l_dic = [0, 1, 2, 3, 4]
eps_dic = []
problem = QuadraticProgramming(G_mat, c_vec, A_mat, b_vec, eps_dic, l_dic)
W_0 = [2]
x0 = np.array([2, 0])
res = active_set(problem, x0, W_0, trajectory=True, tol=1e-9)
xf, Wf, x_traj, q_traj = res

-----------------------------------------------------------
ITERATION 0
-----------------------------------------------------------

Branch 1; k = 0
W0 = [2]
null space used
norm(d0) = 0.2
alpha tilde = 10.0
x1 = [2.2 0.1]
q1 = -0.05
W1 = W0

Branch 2
m2 = -2.3999999999999995 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 1
-----------------------------------------------------------

Branch 1; k = 0
W0 = []
non restriction optimization
norm(d0) = 2.4
alpha tilde = 0.66666667
x1 = [1.4 1.7]
q1 = -6.45
W1 = W0 U [0]

Branch 1; k = 1
W1 = [0]
null space used
norm(d1) = 0.0

Branch 2
mj >= 0 forall j. Solution found.

-----------------------------------------------------------
ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND.
-----------------------------------------------------------

x* = [1.4 1.7]
W* = [0]

Resultados

print("Optimum reached at x* =",xf)
print("Function at x* is q(x*) =",problem.evaluate(xf))
print("Final active set is W* =",Wf)

Optimum reached at x* = [1.4 1.7]
Function at x* is q(x*) = -6.45
Final active set is W* = [0]

x = np.linspace(0, 6, 100)
y = np.linspace(0, 6, 100)
x, y = np.meshgrid(x, y)
z = 0.5 * (problem.G[0,0]*x**2 + 2*problem.G[0,1]*x*y + problem.G[1,1]*y**2) + problem.c[0]*x + problem.c[1]*y

points_x = x_traj[:,0]
points_y = x_traj[:,1]
points_z = q_traj

fig = go.Figure()
fig.add_trace(go.Surface(z=z, x=x, y=y, colorscale='Viridis', opacity=0.9))

fig.add_trace(go.Scatter3d(
    x=points_x,
    y=points_y,
    z=points_z,
    mode='markers+text',
    marker=dict(size=6, color='red'),
    text=[f"P{i}" for i in range (points_x.size)],
    textposition="top center",
    textfont=dict(color="red")
))

fig.update_layout(
    title='Trajectory plot',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z'
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)

fig.show()

A = problem.A
b = problem.b

xtraj = x_traj[:, 0]
ytraj = x_traj[:, 1]

x = np.linspace(0, 6, 900)
y = np.linspace(0, 6, 900)
X, Y = np.meshgrid(x, y)
points = np.vstack((X.ravel(), Y.ravel())).T

feasible = np.all(A @ points.T <= b[:, np.newaxis], axis=0)
X_feasible = X.ravel()[feasible]
Y_feasible = Y.ravel()[feasible]

plt.figure(figsize=(8, 8))
plt.scatter(X_feasible, Y_feasible, color='grey', s=0.5, alpha=0.075)
plt.scatter(xtraj, ytraj, color='red', s=30)
for i in range (len(xtraj)):
    plt.text(xtraj[i], ytraj[i], f"P{i}", fontsize=12, ha='right')

plt.xlim(0, 5)
plt.ylim(0, 2.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Feasible Region and Point Trajectory')
plt.grid(True)
plt.show()

Usamos la biblioteca CVXPY para comparar nuestros resultados.

# Compare the CVXPY solution.
n = int(problem.G.shape[1])
x = cp.Variable(n)
objective = cp.Minimize((1 / 2) * cp.quad_form(x, problem.G) + problem.c.T @ x)
constraints = [problem.A[problem.restrictions_l, :] @ x <= problem.b[problem.restrictions_l]]
prob = cp.Problem(objective, constraints)
prob.solve()
solution = x.value

print("Method used:",prob.solver_stats.solver_name)
print("Our optimum value:", round(float(problem.evaluate(xf)),2))
print("Optimum value with cvxpy:",round(problem.evaluate(x.value),2))
print("Solution x with active set method:",np.round(xf,2))
print("Solution x with cvxpy: ",np.round(solution,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-x.value)/(np.linalg.norm(xf)+np.linalg.norm(x.value)),2))

Method used: OSQP
Our optimum value: -6.45
Optimum value with cvxpy: -6.45
Solution x with active set method: [1.4 1.7]
Solution x with cvxpy:  [1.4 1.7]
Normalized euclidean distance: 0.0

Observaciones

• Las soluciones obtenidas con CVXPY y nuestro implementación del Método del Conjunto Activo coinciden.
• El Método de Conjunto Activo encuentra nuestra solución en dos iteraciones (trés puntos en total).
• CVXPY utiliza el método OSQP, basado en separación de operadores utilizando ADMM (Alternating Direction Method of Multipliers).

Segunda prueba de optimización

Partimos del siguiente punto inicial:

$$ \vec{x}_0 = \begin{pmatrix} 3 \\ 1 \end{pmatrix} $$

y del conjunto activo inicial:

$$ W_0 = \emptyset\subset \{1, 2, 3, 4, 5\} $$

G_mat = np.array([[2, 0], [0, 2]])
A_mat = np.array([[-1, 2], [1, 2], [1, -2], [-1, 0], [0, -1]])
c_vec = np.array([-2, -5])
b_vec = np.array([2, 6, 2, 0, 0])
l_dic = [0, 1, 2, 3, 4]
eps_dic = []
problem = QuadraticProgramming(G_mat, c_vec, A_mat, b_vec, eps_dic, l_dic)
W_0 = [1]
x0 = np.array([3, 1])
res = active_set(problem, x0, W_0, trajectory=True, tol=1e-9)
xf, Wf, x_traj, q_traj = res

-----------------------------------------------------------
ITERATION 0
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1]
null space used
norm(d0) = 2.2
alpha tilde = 0.68181818
x1 = [1.5  1.75]
q1 = -6.437
W1 = W0 U [0]

Branch 1; k = 1
W1 = [1, 0]
null space used
norm(d1) = 0.0

Branch 2
m1 = -0.1250000000000004 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 1
-----------------------------------------------------------

Branch 1; k = 0
W0 = [0]
null space used
norm(d0) = 0.1
alpha tilde = 15.0
x1 = [1.4 1.7]
q1 = -6.45
W1 = W0

Branch 2
mj >= 0 forall j. Solution found.

-----------------------------------------------------------
ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND.
-----------------------------------------------------------

x* = [1.4 1.7]
W* = [0]

Resultados

print("Optimum reached at x* =",xf)
print("Function at x* is q(x*) =",round(problem.evaluate(xf),2))
print("Final active set is W* =",Wf)

Optimum reached at x* = [1.4 1.7]
Function at x* is q(x*) = -6.45
Final active set is W* = [0]

x = np.linspace(0, 6, 100)
y = np.linspace(0, 6, 100)
x, y = np.meshgrid(x, y)
z = 0.5 * (problem.G[0,0]*x**2 + 2*problem.G[0,1]*x*y + problem.G[1,1]*y**2) + problem.c[0]*x + problem.c[1]*y

points_x = x_traj[:,0]
points_y = x_traj[:,1]
points_z = q_traj

fig = go.Figure()
fig.add_trace(go.Surface(z=z, x=x, y=y, colorscale='Viridis', opacity=0.9))

fig.add_trace(go.Scatter3d(
    x=points_x,
    y=points_y,
    z=points_z,
    mode='markers+text',
    marker=dict(size=6, color='red'),
    text=[f"P{i}" for i in range (len(points_x))],
    textposition="top center",
    textfont=dict(color="red")
))

fig.update_layout(
    title='Trajectory plot',
    scene=dict(
        xaxis_title='X',
        yaxis_title='Y',
        zaxis_title='Z'
    ),
    margin=dict(l=0, r=0, b=0, t=30)
)

fig.show()

A = problem.A
b = problem.b

xtraj = x_traj[:, 0]
ytraj = x_traj[:, 1]

x = np.linspace(0, 6, 900)
y = np.linspace(0, 6, 900)
X, Y = np.meshgrid(x, y)
points = np.vstack((X.ravel(), Y.ravel())).T

feasible = np.all(A @ points.T <= b[:, np.newaxis], axis=0)
X_feasible = X.ravel()[feasible]
Y_feasible = Y.ravel()[feasible]

plt.figure(figsize=(8, 8))
plt.scatter(X_feasible, Y_feasible, color='grey', s=0.5, alpha=0.075)
plt.scatter(xtraj, ytraj, color='red', s=30)
for i in range (len(xtraj)):
    plt.text(xtraj[i], ytraj[i], f"P{i}", fontsize=12, ha='right')

plt.xlim(0, 5)
plt.ylim(0, 2.5)
plt.xlabel('x')
plt.ylabel('y')
plt.title('Feasible Region and Point Trajectory')
plt.grid(True)
plt.show()

# Compare with CVXPY solution
n = int(problem.G.shape[1])
x = cp.Variable(n)
objective = cp.Minimize((1 / 2) * cp.quad_form(x, problem.G) + problem.c.T @ x)
constraints = [problem.A[problem.restrictions_l, :] @ x <= problem.b[problem.restrictions_l]]
prob = cp.Problem(objective, constraints)
prob.solve()
solution = x.value

print("Our optimum value:", round(float(problem.evaluate(xf)),2))
print("Optimum value obtained with cvxpy:",round(problem.evaluate(x.value),2))
print("Solution x with active set method:",np.round(xf,2))
print("Solution x with cvxpy: ",np.round(solution,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-x.value)/(np.linalg.norm(xf)+np.linalg.norm(x.value)),2))

Our optimum value: -6.45
Optimum value obtained with cvxpy: -6.45
Solution x with active set method: [1.4 1.7]
Solution x with cvxpy:  [1.4 1.7]
Normalized euclidean distance: 0.0

Observaciones

• Obtenemos una solución idéntica
• En este caso, comenzamos en el interior de nuestra región factible. El algoritmo se mueve a la frontera de la región y finalmente optimiza en la frontera.

---

IV. Problema - Klee Minty

Para este problema, se considera la matriz $ G = 10^{-4} I $ y el problema cuadrático de Klee-Minty definido por:

$$ \min_{\vec{x} \in \mathbb{R}^n} \left\{ \frac{1}{2} \vec{x}^\top G \vec{x} - \sum_{i=1}^n x_i \right\} $$

sujeto a las siguientes restricciones:

$$ \begin{aligned} x_1 &\leq 1 \\ 2 \sum_{j=1}^{i-1} x_j + x_i &\leq 2^{i-1}, \quad i = 2, \dots, n \\ x_1, \dots, x_n &\geq 0 \end{aligned} $$

Sea $n = 20$. Este problema impone $n$ restricciones de no negatividad, una para cada variable. Para aplicar el método de optimización (por ejemplo, un método de conjunto activo o de proyecciones), se debe construir un punto inicial $\vec{x}_{0}$ que cumpla con un conjunto activo inicial $W_0$. Este conjunto debe ser un subconjunto aleatorio formado por 8 de las últimas 10 restricciones de no negatividad.

A partir del punto $\vec{x}_{0}$ y del conjunto activo $W_0$, se debe ejecutar el método seleccionado y determinar:

• El conjunto activo inicial $W_0$,
• El número de iteraciones necesarias,
• El valor óptimo alcanzado.

def klee_minty(n):
    G = np.identity(n) * 1e-4
    c = np.ones(n) * -1
    b = 2 ** np.arange(1, n + 1) - 1
    b = np.concatenate((b, np.zeros(n)))
    L = np.tril(np.ones((n, n))) + np.tril(np.ones((n, n)), k=-1)
    Identity = np.identity(n)
    A = np.vstack((L, -Identity))
    restrictions_eps = []
    restrictions_l = [i for i in range(n)]
    problem = QuadraticProgramming(G, c, A, b, restrictions_eps, restrictions_l)
    return problem


problem_km = klee_minty(20)
k = random.randint(1, 8)
W0 = (
    np.random.choice(problem_km.restrictions_l[-8:], k, replace=False)
    .astype(int)
    .tolist()
)
problem_km_0=klee_minty(20)
problem_km_0.restrictions_eps = list(set(problem_km_0.restrictions_eps + W0))
problem_km_0.restrictions_l = [
    i for i in problem_km_0.restrictions_l if i not in W0
]
x0 = problem_km_0.generate_feasible()
res = active_set(problem_km, x0, W0, max_iter=100000,
                 tol=1e-9, trajectory=True)
xf, Wf, x_traj, q_traj = res

-----------------------------------------------------------
ITERATION 0
-----------------------------------------------------------

Branch 1; k = 0
W0 = [18, 12, 17, 14, 13, 16, 19, 15]
schur complement used
norm(d0) = 23892.447
No alpha tilde.
x1 = [  -952.852   -952.852   -952.852   -952.852   -952.852   -952.852
   -952.852   -952.852   -952.852   -952.852   -952.852   -952.852
  31059.447 -22867.447  39251.447  -6483.447  72019.447  59052.553
 203091.447 321196.553]
q1 = 7123543.321
W1 = W0

Branch 2
m19 = -31.119655324675243 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 1
-----------------------------------------------------------

Branch 1; k = 0
W0 = [18, 12, 17, 14, 13, 16, 15]
schur complement used
norm(d0) = 311196.553
alpha tilde = 0.51647119
x1 = [1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00
 1.000000e+00 1.000000e+00 8.167000e+03 2.500000e+01 1.635900e+04
 1.640900e+04 4.912700e+04 8.194500e+04 1.801990e+05 1.604725e+05]
q1 = 2885040.535
W1 = W0 U [1]

Branch 1; k = 1
W1 = [18, 12, 17, 14, 13, 16, 15, 1]
schur complement used
norm(d1) = 150472.5
alpha tilde = 0.00136048
x2 = [7.14000000e-01 1.57100000e+00 2.42900000e+00 2.42900000e+00
 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00
 2.42900000e+00 2.42900000e+00 2.42900000e+00 2.42900000e+00
 8.13785700e+03 5.41430000e+01 1.63298570e+04 1.64381430e+04
 4.90978570e+04 8.19741430e+04 1.80169857e+05 1.60267786e+05]
q2 = 2881524.022
W2 = W1 U [2]

Branch 1; k = 2
W2 = [18, 12, 17, 14, 13, 16, 15, 1, 2]
schur complement used
norm(d2) = 150267.786
alpha tilde = 0.00219133
x3 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 5.000000e+00
 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00 5.000000e+00
 5.000000e+00 5.000000e+00 8.091000e+03 1.010000e+02 1.628300e+04
 1.648500e+04 4.905100e+04 8.202100e+04 1.801230e+05 1.599385e+05]
q3 = 2875877.999
W3 = W2 U [3]

Branch 1; k = 3
W3 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3]
null space used
norm(d3) = 149938.5
alpha tilde = 0.00393895
x4 = [6.000000e-01 1.800000e+00 2.200000e+00 5.800000e+00 1.020000e+01
 1.020000e+01 1.020000e+01 1.020000e+01 1.020000e+01 1.020000e+01
 1.020000e+01 1.020000e+01 8.007000e+03 1.850000e+02 1.619900e+04
 1.656900e+04 4.896700e+04 8.210500e+04 1.800390e+05 1.593479e+05]
q4 = 2865783.069
W4 = W3 U [4]

Branch 1; k = 4
W4 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4]
null space used
norm(d4) = 149347.9
alpha tilde = 0.00452233
x5 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 1.700000e+01 1.700000e+01 1.700000e+01 1.700000e+01 1.700000e+01
 1.700000e+01 1.700000e+01 7.911000e+03 2.810000e+02 1.610300e+04
 1.666500e+04 4.887100e+04 8.220100e+04 1.799430e+05 1.586725e+05]
q5 = 2854288.422
W5 = W4 U [0]

Branch 1; k = 5
W5 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0]
null space used
norm(d5) = 148672.5
alpha tilde = 0.00265012
x6 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 2.100000e+01 2.100000e+01 2.100000e+01 2.100000e+01
 2.100000e+01 2.100000e+01 7.855000e+03 3.370000e+02 1.604700e+04
 1.672100e+04 4.881500e+04 8.225700e+04 1.798870e+05 1.582785e+05]
q6 = 2847607.025
W6 = W5 U [5]

Branch 1; k = 6
W6 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5]
null space used
norm(d6) = 148278.5
alpha tilde = 0.01253722
x7 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 4.300000e+01 4.300000e+01 4.300000e+01 4.300000e+01
 4.300000e+01 4.300000e+01 7.591000e+03 6.010000e+02 1.578300e+04
 1.698500e+04 4.855100e+04 8.252100e+04 1.796230e+05 1.564195e+05]
q7 = 2816324.832
W7 = W6 U [6]

Branch 1; k = 7
W7 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6]
null space used
norm(d7) = 146419.5
alpha tilde = 0.02022272
x8 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 4.300000e+01 8.500000e+01 8.500000e+01 8.500000e+01
 8.500000e+01 8.500000e+01 7.171000e+03 1.021000e+03 1.536300e+04
 1.740500e+04 4.813100e+04 8.294100e+04 1.792030e+05 1.534585e+05]
q8 = 2767320.956
W8 = W7 U [7]

Branch 1; k = 8
W8 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7]
null space used
norm(d8) = 143458.5
alpha tilde = 0.03387042
x9 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 1.710000e+02
 1.710000e+02 1.710000e+02 6.483000e+03 1.709000e+03 1.467500e+04
 1.809300e+04 4.744300e+04 8.362900e+04 1.785150e+05 1.485995e+05]
q9 = 2689092.17
W9 = W8 U [8]

Branch 1; k = 9
W9 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8]
null space used
norm(d9) = 138599.5
alpha tilde = 0.05212862
x10 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02
 3.410000e+02 3.410000e+02 5.463000e+03 2.729000e+03 1.365500e+04
 1.911300e+04 4.642300e+04 8.464900e+04 1.774950e+05 1.413745e+05]
q10 = 2577795.388
W10 = W9 U [9]

Branch 1; k = 10
W10 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8, 9]
null space used
norm(d10) = 131374.5
alpha tilde = 0.07419248
x11 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02
 6.830000e+02 6.830000e+02 4.095000e+03 4.097000e+03 1.228700e+04
 2.048100e+04 4.505500e+04 8.601700e+04 1.761270e+05 1.316275e+05]
q11 = 2437189.527
W11 = W10 U [10]

Branch 1; k = 11
W11 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8, 9, 10]
null space used
norm(d11) = 121627.5
alpha tilde = 0.08130563
x12 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02
 6.830000e+02 1.365000e+03 2.731000e+03 5.461000e+03 1.092300e+04
 2.184500e+04 4.369100e+04 8.738100e+04 1.747630e+05 1.217385e+05]
q12 = 2305886.147
W12 = W11 U [11]

Branch 1; k = 12
W12 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 0, 5, 6, 7, 8, 9, 10, 11]
null space used
norm(d12) = 111738.5
No alpha tilde.
x13 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01
 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03
 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 4.36910e+04 8.73810e+04
 1.74763e+05 1.00000e+04]
q13 = 1681611.528
W13 = W12

Branch 2
m0 = -22.302900000007853 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 2
-----------------------------------------------------------

Branch 1; k = 0
W0 = [18, 12, 17, 14, 13, 16, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
null space used
norm(d0) = 6110.384
No alpha tilde.
x1 = [ -3054.192   6111.384  -6107.384   6115.384  -6099.384   6131.384
  -6067.384   6195.384  -5939.384   6451.384  -5427.384   7475.384
  -3379.384  11571.384   4812.616  27955.384  37580.616  93491.384
 168652.616  10000.   ]
q1 = 1647541.709
W1 = W0

Branch 2
m16 = -17.79030821918008 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 3
-----------------------------------------------------------

Branch 1; k = 0
W0 = [18, 12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11]
null space used
norm(d0) = 40243.785
alpha tilde = 0.95655556
x1 = [1.000000e+00 1.000000e+00 3.000000e+00 5.000000e+00 1.100000e+01
 2.100000e+01 4.300000e+01 8.500000e+01 1.710000e+02 3.410000e+02
 6.830000e+02 1.365000e+03 2.731000e+03 5.461000e+03 1.092300e+04
 2.184500e+04 2.138810e+04 1.319868e+05 1.301572e+05 1.000000e+04]
q1 = 1440535.133
W1 = W0 U [0]

Branch 1; k = 1
W1 = [18, 12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0]
null space used
norm(d1) = 1717.533
No alpha tilde.
x2 = [1.00000000e+00 1.00000000e+00 3.00000000e+00 5.00000000e+00
 1.10000000e+01 2.10000000e+01 4.30000000e+01 8.50000000e+01
 1.71000000e+02 3.41000000e+02 6.83000000e+02 1.36500000e+03
 2.73100000e+03 5.46100000e+03 1.09230000e+04 2.18450000e+04
 2.05293330e+04 1.33704333e+05 1.28439667e+05 1.00000000e+04]
q2 = 1440203.266
W2 = W1

Branch 2
m18 = -11.843966666666557 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 4
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0]
null space used
norm(d0) = 118439.667
alpha tilde = 0.48889167
x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01
 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03
 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 4.36910e+04 8.73810e+04
 7.05355e+04 1.00000e+04]
q1 = 507496.557
W1 = W0 U [16]

Branch 1; k = 1
W1 = [12, 17, 14, 13, 15, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16]
null space used
norm(d1) = 60535.5
No alpha tilde.
x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01
 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03
 2.7310e+03 5.4610e+03 1.0923e+04 2.1845e+04 4.3691e+04 8.7381e+04
 1.0000e+04 1.0000e+04]
q2 = 324269.219
W2 = W1

Branch 2
m3 = -10.649699999991908 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 5
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 17, 14, 13, 15, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16]
null space used
norm(d0) = 3736.737
No alpha tilde.
x1 = [ 1.0000000e+00  1.0000000e+00  3.0000000e+00 -1.8633680e+03
  3.7477370e+03 -3.7157370e+03  3.7797370e+03 -3.6517370e+03
  3.9077370e+03 -3.3957370e+03  4.4197370e+03 -2.3717370e+03
  6.4677370e+03  1.7242630e+03  1.4659737e+04  1.8108263e+04
  4.7427737e+04  8.3644263e+04  1.0000000e+04  1.0000000e+04]
q1 = 314320.438
W1 = W0

Branch 2
m15 = -8.054131578943784 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 6
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 17, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16]
null space used
norm(d0) = 18331.437
alpha tilde = 0.95806252
x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01
 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03
 2.73100e+03 5.46100e+03 1.09230e+04 1.11953e+04 6.49904e+04 6.60816e+04
 1.00000e+04 1.00000e+04]
q1 = 269634.821
W1 = W0 U [3]

Branch 1; k = 1
W1 = [12, 17, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16, 3]
null space used
norm(d1) = 750.6
No alpha tilde.
x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01
 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03
 2.7310e+03 5.4610e+03 1.0923e+04 1.0820e+04 6.5741e+04 6.5331e+04
 1.0000e+04 1.0000e+04]
q2 = 269571.438
W2 = W1

Branch 2
m17 = -5.533100000000015 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 7
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16, 3]
null space used
norm(d0) = 55331.0
alpha tilde = 0.49813848
x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01
 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03
 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 4.36910e+04 3.77685e+04
 1.00000e+04 1.00000e+04]
q1 = 63432.741
W1 = W0 U [15]

Branch 1; k = 1
W1 = [12, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 16, 3, 15]
null space used
norm(d1) = 27768.5
No alpha tilde.
x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01
 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03
 2.7310e+03 5.4610e+03 1.0923e+04 2.1845e+04 4.3691e+04 1.0000e+04
 1.0000e+04 1.0000e+04]
q2 = 24878.261
W2 = W1

Branch 2
m0 = -4.826500000005277 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 8
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 14, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 16, 3, 15]
null space used
norm(d0) = 1485.077
No alpha tilde.
x1 = [ -741.538  1486.077 -1482.077  1490.077 -1474.077  1506.077 -1442.077
  1570.077 -1314.077  1826.077  -802.077  2850.077  1245.923  6946.077
  9437.923 23330.077 42205.923 10000.    10000.    10000.   ]
q1 = 23086.33
W1 = W0

Branch 2
m14 = -3.7189615384644465 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 9
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 16, 3, 15]
null space used
norm(d0) = 8434.758
alpha tilde = 0.96836482
x1 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01
 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03
 2.7310e+03 5.4610e+03 6.0965e+03 3.1498e+04 3.4038e+04 1.0000e+04
 1.0000e+04 1.0000e+04]
q1 = 13827.628
W1 = W0 U [0]

Branch 1; k = 1
W1 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 16, 3, 15, 0]
null space used
norm(d1) = 261.444
No alpha tilde.
x2 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01
 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02
 6.8300000e+02 1.3650000e+03 2.7310000e+03 5.4610000e+03 5.9657780e+03
 3.1759444e+04 3.3776556e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04]
q2 = 13819.938
W2 = W1

Branch 2
m16 = -2.3776555555555046 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 10
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 3, 15, 0]
null space used
norm(d0) = 23776.556
alpha tilde = 0.52123006
x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01
 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03
 2.73100e+03 5.46100e+03 1.09230e+04 2.18450e+04 2.13835e+04 1.00000e+04
 1.00000e+04 1.00000e+04]
q1 = -25396.709
W1 = W0 U [14]

Branch 1; k = 1
W1 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 3, 15, 0, 14]
null space used
norm(d1) = 11383.5
No alpha tilde.
x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01
 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03
 2.7310e+03 5.4610e+03 1.0923e+04 2.1845e+04 1.0000e+04 1.0000e+04
 1.0000e+04 1.0000e+04]
q2 = -31875.913
W2 = W1

Branch 2
m3 = -1.9116999999972746 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 11
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 13, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 15, 0, 14]
null space used
norm(d0) = 780.286
No alpha tilde.
x1 = [ 1.0000000e+00  1.0000000e+00  3.0000000e+00 -3.8514300e+02
  7.9128600e+02 -7.5928600e+02  8.2328600e+02 -6.9528600e+02
  9.5128600e+02 -4.3928600e+02  1.4632860e+03  5.8471400e+02
  3.5112860e+03  4.6807140e+03  1.1703286e+04  2.1064714e+04
  1.0000000e+04  1.0000000e+04  1.0000000e+04  1.0000000e+04]
q1 = -32248.831
W1 = W0

Branch 2
m13 = -1.3403571428555356 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 12
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 15, 0, 14]
null space used
norm(d0) = 3065.92
alpha tilde = 0.99256168
x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01
 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03
 2.73100e+03 3.54930e+03 1.47464e+04 1.80216e+04 1.00000e+04 1.00000e+04
 1.00000e+04 1.00000e+04]
q1 = -33539.541
W1 = W0 U [3]

Branch 1; k = 1
W1 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 15, 0, 14, 3]
null space used
norm(d1) = 22.156
No alpha tilde.
x2 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01
 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02
 6.8300000e+02 1.3650000e+03 2.7310000e+03 3.5382220e+03 1.4768556e+04
 1.7999444e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04]
q2 = -33539.596
W2 = W1

Branch 2
m15 = -0.7999444444444724 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 13
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 14, 3]
null space used
norm(d0) = 7999.444
alpha tilde = 0.60090979
x1 = [1.00000e+00 1.00000e+00 3.00000e+00 5.00000e+00 1.10000e+01 2.10000e+01
 4.30000e+01 8.50000e+01 1.71000e+02 3.41000e+02 6.83000e+02 1.36500e+03
 2.73100e+03 5.46100e+03 1.09230e+04 1.31925e+04 1.00000e+04 1.00000e+04
 1.00000e+04 1.00000e+04]
q1 = -38381.511
W1 = W0 U [13]

Branch 1; k = 1
W1 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 0, 14, 3, 13]
null space used
norm(d1) = 3192.5
No alpha tilde.
x2 = [1.0000e+00 1.0000e+00 3.0000e+00 5.0000e+00 1.1000e+01 2.1000e+01
 4.3000e+01 8.5000e+01 1.7100e+02 3.4100e+02 6.8300e+02 1.3650e+03
 2.7310e+03 5.4610e+03 1.0923e+04 1.0000e+04 1.0000e+04 1.0000e+04
 1.0000e+04 1.0000e+04]
q2 = -38891.114
W2 = W1

Branch 2
m0 = -0.45730000000121807 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 14
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 10, 11, 14, 3, 13]
null space used
norm(d0) = 160.456
No alpha tilde.
x1 = [-7.9228000e+01  1.6145600e+02 -1.5745600e+02  1.6545600e+02
 -1.4945600e+02  1.8145600e+02 -1.1745600e+02  2.4545600e+02
  1.0544000e+01  5.0145600e+02  5.2254400e+02  1.5254560e+03
  2.5705440e+03  5.6214560e+03  1.0762544e+04  1.0000000e+04
  1.0000000e+04  1.0000000e+04  1.0000000e+04  1.0000000e+04]
q1 = -38909.458
W1 = W0

Branch 2
m10 = -0.2894894736850741 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 15
-----------------------------------------------------------

Branch 1; k = 0
W0 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13]
null space used
norm(d0) = 350.08
alpha tilde = 0.9930721
x1 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01
 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02
 4.2894400e+02 1.8731110e+03 2.2228890e+03 5.9691110e+03 1.0414889e+04
 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04]
q1 = -38946.486
W1 = W0 U [0]

Branch 1; k = 1
W1 = [12, 1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13, 0]
null space used
norm(d1) = 2.359
No alpha tilde.
x2 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01
 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02
 4.2776500e+02 1.8754710e+03 2.2205290e+03 5.9714710e+03 1.0412529e+04
 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04]
q2 = -38946.487
W2 = W1

Branch 2
m12 = -0.11026470588274505 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 16
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13, 0]
null space used
norm(d0) = 416.093
alpha tilde = 1.22681782
x1 = [1.0000000e+00 1.0000000e+00 3.0000000e+00 5.0000000e+00 1.1000000e+01
 2.1000000e+01 4.3000000e+01 8.5000000e+01 1.7100000e+02 3.4100000e+02
 6.3581100e+02 1.4593770e+03 2.2829430e+03 6.2627360e+03 1.0121264e+04
 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04 1.0000000e+04]
q1 = -38965.987
W1 = W0

Branch 2
m0 = -0.018681132075438582 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 17
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1, 2, 4, 5, 6, 7, 8, 9, 11, 14, 3, 13]
null space used
norm(d0) = 9.876
alpha tilde = 8.16878223
x1 = [-3.9380000e+00  1.0876000e+01 -6.8760000e+00  1.4876000e+01
  1.1240000e+00  3.0876000e+01  3.3124000e+01  9.4876000e+01
  1.6112400e+02  3.5087600e+02  6.3171200e+02  1.4577000e+03
  2.2836890e+03  6.2629220e+03  1.0121078e+04  1.0000000e+04
  1.0000000e+04  1.0000000e+04  1.0000000e+04  1.0000000e+04]
q1 = -38966.033
W1 = W0

Branch 2
m14 = -0.012107780548628416 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 18
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1, 2, 4, 5, 6, 7, 8, 9, 11, 3, 13]
null space used
norm(d0) = 121.078
alpha tilde = 6.06203426
x1 = [-3.717000e+00  1.043500e+01 -6.435000e+00  1.443500e+01  1.565000e+00
  3.043500e+01  3.356500e+01  9.443500e+01  1.615650e+02  3.504350e+02
  6.235440e+02  1.474477e+03  2.325409e+03  6.162705e+03  1.000000e+04
  1.000000e+04  1.000000e+04  1.000000e+04  1.000000e+04  1.000000e+04]
q1 = -38967.372
W1 = W0

Branch 2
m6 = -0.00905004557884357 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 19
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1, 2, 4, 5, 7, 8, 9, 11, 3, 13]
null space used
norm(d0) = 13.73
alpha tilde = 1.06591256
x1 = [7.080000e-01 1.583000e+00 2.417000e+00 5.583000e+00 1.041700e+01
 2.158300e+01 3.112600e+01 1.081650e+02 1.478350e+02 3.641650e+02
 6.178630e+02 1.472110e+03 2.326356e+03 6.163178e+03 1.000000e+04
 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04]
q1 = -38967.423
W1 = W0

Branch 2
m8 = -0.0037368474923263548 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 20
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1, 2, 4, 5, 7, 9, 11, 3, 13]
schur complement used
norm(d0) = 15.423
alpha tilde = 1.42652892
x1 = [3.64000e-01 2.27200e+00 1.72800e+00 6.27200e+00 9.72800e+00 2.22720e+01
 3.83520e+01 9.30230e+01 1.47694e+02 3.79588e+02 6.11481e+02 1.46945e+03
 2.32742e+03 6.16371e+03 1.00000e+04 1.00000e+04 1.00000e+04 1.00000e+04
 1.00000e+04 1.00000e+04]
q1 = -38967.452
W1 = W0

Branch 2
m4 = -0.0003537463373766392 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 21
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1, 2, 5, 7, 9, 11, 3, 13]
schur complement used
norm(d0) = 1.317
alpha tilde = 2.12154936
x1 = [6.640000e-01 1.672000e+00 2.328000e+00 5.672000e+00 9.370000e+00
 2.358800e+01 3.780700e+01 9.279800e+01 1.477880e+02 3.796260e+02
 6.114650e+02 1.469444e+03 2.327423e+03 6.163711e+03 1.000000e+04
 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04]
q1 = -38967.452
W1 = W0

Branch 2
m2 = -3.5290783356346474e-05 < 0. Back to Branch 1.

-----------------------------------------------------------
ITERATION 22
-----------------------------------------------------------

Branch 1; k = 0
W0 = [1, 5, 7, 9, 11, 3, 13]
schur complement used
norm(d0) = 0.145
alpha tilde = 4.73835548
x1 = [7.350000e-01 1.530000e+00 2.326000e+00 5.818000e+00 9.309000e+00
 2.356300e+01 3.781700e+01 9.280200e+01 1.477860e+02 3.796260e+02
 6.114650e+02 1.469444e+03 2.327422e+03 6.163711e+03 1.000000e+04
 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04]
q1 = -38967.452
W1 = W0

Branch 2
mj >= 0 forall j. Solution found.

-----------------------------------------------------------
ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND.
-----------------------------------------------------------

x* = [7.350000e-01 1.530000e+00 2.326000e+00 5.818000e+00 9.309000e+00
 2.356300e+01 3.781700e+01 9.280200e+01 1.477860e+02 3.796260e+02
 6.114650e+02 1.469444e+03 2.327422e+03 6.163711e+03 1.000000e+04
 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04]
W* = [1, 5, 7, 9, 11, 3, 13]

Resultados

print("Optimum reached at x* =",np.round(xf,decimals=3))
print("Function at x* is q(x*) =",np.round(problem_km.evaluate(xf),decimals=3))
print("Final active set is W* =",Wf)

Optimum reached at x* = [7.350000e-01 1.530000e+00 2.326000e+00 5.818000e+00 9.309000e+00
 2.356300e+01 3.781700e+01 9.280200e+01 1.477860e+02 3.796260e+02
 6.114650e+02 1.469444e+03 2.327422e+03 6.163711e+03 1.000000e+04
 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04 1.000000e+04]
Function at x* is q(x*) = -38967.452
Final active set is W* = [1, 5, 7, 9, 11, 3, 13]

n = q_traj.size
x = np.arange(n)

fig, axs = plt.subplots(1, 1, figsize=(12, 5))
fig.suptitle('q optimization per iteration', fontsize=14)
axs.plot(x, q_traj, marker='o', color='black')
axs.set_xlabel('Step')
axs.set_xticks(x)
axs.set_ylabel('q')
axs.grid(True)

n = int(problem_km.G.shape[1])
x = cp.Variable(n)
prob = cp.Problem(
    cp.Minimize((1 / 2) * cp.quad_form(x, problem_km.G) +
                problem_km.c.T @ x), [problem_km.A[problem_km.restrictions_l,:] @ x <= problem_km.b[problem_km.restrictions_l]]
)
prob.solve()
solution_cvxpy = x.value

print("Method used:",prob.solver_stats.solver_name)
print("Our optimum value:", round(float(problem_km.evaluate(xf)),2))
print("Optimum value obtained with cvxpy:",round(prob.value,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-solution_cvxpy)/(np.linalg.norm(xf)+np.linalg.norm(solution_cvxpy)),2))

Method used: OSQP
Our optimum value: -38967.45
Optimum value obtained with cvxpy: -38967.45
Normalized euclidean distance: 0.0

Observaciones

• La solución obtenida con nuestro algoritmo es identica a la obtenida con CVXPY (distancia normalizada cerca de 0).
• Utilizamos x0 inicial factible e imponiendo a W0 como restricciones activas.

print("W0:",W0)
print("iteraciones:",x_traj.size)
print("Valor óptimo obtenido:",round(float(problem_km.evaluate(xf)),2))

W0: [18, 12, 17, 14, 13, 16, 19, 15]
iteraciones: 1340
Valor óptimo obtenido: -38967.45

---

V. Problema 3 - lp_blend

A partir de los datos del archivo lp_blend.mat, se plantea el siguiente problema de optimización cuadrática. Dado $G \in \mathbb{R}^{n \times n}$, se desea resolver:

$$ \begin{aligned} \text{minimizar} \quad & \vec{x}^\top G \vec{x} + \vec{c}^\top \vec{x} \\ \text{sujeto a} \quad & A\vec{x} = \vec{b}, \\ & \vec{x} \geq 0. \end{aligned} $$

Este planteamiento puede implementarse dentro de una función separada. Una vez definido el problema, se debe aplicar el algoritmo desarrollado para resolverlo, siguiendo los pasos:

• Encontrar un punto inicial $\vec{x_0}$ mediante programación lineal (utilizando from scipy.optimize import linprog), considerando solo $A$ y $\vec{b}$.
• Definir el conjunto activo inicial $W_0$ y un vector $z \in E$.

Se deben documentar los siguientes elementos:

• El conjunto activo inicial $W_0$,
• El número de iteraciones realizadas,
• El valor óptimo alcanzado.

problem_load = loadProblem("lp_blend.mat")
Ae = problem_load["AE"]
n = Ae.shape[1]
Al = -np.identity(n)
A = np.vstack((Ae, Al))
be = problem_load["bE"]
bl = np.zeros(n)
b = np.concatenate((be, bl))
c = problem_load["c"]
G = np.identity(n)
restrictions_eps = np.arange(0, Ae.shape[0]).tolist()
restrictions_l = np.arange(Ae.shape[0], Ae.shape[0] + n).tolist()
problem_blend = QuadraticProgramming(G, c, A, b, restrictions_eps, restrictions_l)
x0 = problem_blend.generate_feasible()
W0 = restrictions_eps.copy()
res = active_set(problem_blend, x0, W0, trajectory=True, tol=1e-9)
xf, wf, x_traj, q_traj = res

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf]

-----------------------------------------------------------
ITERATION 0
-----------------------------------------------------------

Branch 1; k = 0
W0 = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73]
null space used
norm(d0) = 14.405
alpha tilde = 1.47629031
x1 = [-3.6200e-01 -1.3760e+00  1.6860e+00  2.3000e-02  1.5700e-01 -4.8600e-01
  5.4000e-02 -3.5400e-01  2.0320e+00  5.1000e-02  1.2100e-01  6.8000e-02
  2.6000e-01  7.3000e-02  1.9570e+00  6.0000e-02 -3.1000e-02  1.8200e-01
  7.6500e-01  1.7830e+00  1.6510e+00  8.8550e+00  1.6940e+00  1.5243e+01
  1.2154e+01  1.0247e+01  3.5510e+00  4.6240e+00  1.1863e+01  1.4000e-01
 -9.8000e-02  1.1077e+01  8.8960e+00  5.9480e+00  2.6070e+00 -3.2500e-01
  2.5000e-02  1.1700e-01  3.1750e+00  5.9000e-02 -6.2800e-01  1.9310e+00
  1.2440e+00  1.9900e+00  6.1600e-01 -6.9500e-01 -3.0300e-01 -1.4420e+00
  1.4690e+00  1.9180e+00  2.3090e+00  7.7600e-01  3.7300e-01 -1.2360e+00
 -6.2700e-01  5.9700e-01  6.7300e-01  4.1400e-01 -1.2000e-01  3.6500e-01
  1.1840e+00  4.6300e-01  1.1810e+00 -1.2360e+00 -6.2700e-01  1.9000e-02
  1.9990e+00  1.2280e+00  4.3700e-01 -2.7100e-01  3.2500e-01  4.1000e-02
  7.4700e-01  6.1000e-02  1.9890e+00  1.3400e+00  8.0500e-01 -4.2500e-01
  1.1200e-01  2.4000e-02  7.0400e-01  1.2000e-02  1.7380e+00  1.2320e+00
  1.7300e-01  1.0450e+00  6.4100e-01  1.6850e+00  4.4700e-01  2.9870e+00
  3.4350e+00  2.5230e+00  2.1190e+00  4.6420e+00  4.1000e-01  1.6010e+00
  1.5500e+00  2.7800e-01  4.1000e-02  1.1700e-01  3.9960e+00  4.2050e+00
 -1.3400e-01 -5.5300e-01 -1.3400e-01  4.6300e-01  1.1600e-01 -2.5000e-01
  2.5200e-01 -2.7830e+00 -4.3000e-01  3.0000e-02  4.6200e+00  1.0070e+00]
q1 = 587.078
W1 = W0

Branch 2
mj >= 0 forall j. Solution found.

-----------------------------------------------------------
ACTIVE SET METHOD SUCCESFUL. SOLUTION FOUND.
-----------------------------------------------------------

x* = [-3.6200e-01 -1.3760e+00  1.6860e+00  2.3000e-02  1.5700e-01 -4.8600e-01
  5.4000e-02 -3.5400e-01  2.0320e+00  5.1000e-02  1.2100e-01  6.8000e-02
  2.6000e-01  7.3000e-02  1.9570e+00  6.0000e-02 -3.1000e-02  1.8200e-01
  7.6500e-01  1.7830e+00  1.6510e+00  8.8550e+00  1.6940e+00  1.5243e+01
  1.2154e+01  1.0247e+01  3.5510e+00  4.6240e+00  1.1863e+01  1.4000e-01
 -9.8000e-02  1.1077e+01  8.8960e+00  5.9480e+00  2.6070e+00 -3.2500e-01
  2.5000e-02  1.1700e-01  3.1750e+00  5.9000e-02 -6.2800e-01  1.9310e+00
  1.2440e+00  1.9900e+00  6.1600e-01 -6.9500e-01 -3.0300e-01 -1.4420e+00
  1.4690e+00  1.9180e+00  2.3090e+00  7.7600e-01  3.7300e-01 -1.2360e+00
 -6.2700e-01  5.9700e-01  6.7300e-01  4.1400e-01 -1.2000e-01  3.6500e-01
  1.1840e+00  4.6300e-01  1.1810e+00 -1.2360e+00 -6.2700e-01  1.9000e-02
  1.9990e+00  1.2280e+00  4.3700e-01 -2.7100e-01  3.2500e-01  4.1000e-02
  7.4700e-01  6.1000e-02  1.9890e+00  1.3400e+00  8.0500e-01 -4.2500e-01
  1.1200e-01  2.4000e-02  7.0400e-01  1.2000e-02  1.7380e+00  1.2320e+00
  1.7300e-01  1.0450e+00  6.4100e-01  1.6850e+00  4.4700e-01  2.9870e+00
  3.4350e+00  2.5230e+00  2.1190e+00  4.6420e+00  4.1000e-01  1.6010e+00
  1.5500e+00  2.7800e-01  4.1000e-02  1.1700e-01  3.9960e+00  4.2050e+00
 -1.3400e-01 -5.5300e-01 -1.3400e-01  4.6300e-01  1.1600e-01 -2.5000e-01
  2.5200e-01 -2.7830e+00 -4.3000e-01  3.0000e-02  4.6200e+00  1.0070e+00]
W* = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73]

Resultados

print("Optimum reached at x* =",np.round(xf,decimals=3))
print("Function at x* is q(x*) =",np.round(problem_blend.evaluate(xf),decimals=3))
print("Final active set is W* =",Wf)

Optimum reached at x* = [-3.6200e-01 -1.3760e+00  1.6860e+00  2.3000e-02  1.5700e-01 -4.8600e-01
  5.4000e-02 -3.5400e-01  2.0320e+00  5.1000e-02  1.2100e-01  6.8000e-02
  2.6000e-01  7.3000e-02  1.9570e+00  6.0000e-02 -3.1000e-02  1.8200e-01
  7.6500e-01  1.7830e+00  1.6510e+00  8.8550e+00  1.6940e+00  1.5243e+01
  1.2154e+01  1.0247e+01  3.5510e+00  4.6240e+00  1.1863e+01  1.4000e-01
 -9.8000e-02  1.1077e+01  8.8960e+00  5.9480e+00  2.6070e+00 -3.2500e-01
  2.5000e-02  1.1700e-01  3.1750e+00  5.9000e-02 -6.2800e-01  1.9310e+00
  1.2440e+00  1.9900e+00  6.1600e-01 -6.9500e-01 -3.0300e-01 -1.4420e+00
  1.4690e+00  1.9180e+00  2.3090e+00  7.7600e-01  3.7300e-01 -1.2360e+00
 -6.2700e-01  5.9700e-01  6.7300e-01  4.1400e-01 -1.2000e-01  3.6500e-01
  1.1840e+00  4.6300e-01  1.1810e+00 -1.2360e+00 -6.2700e-01  1.9000e-02
  1.9990e+00  1.2280e+00  4.3700e-01 -2.7100e-01  3.2500e-01  4.1000e-02
  7.4700e-01  6.1000e-02  1.9890e+00  1.3400e+00  8.0500e-01 -4.2500e-01
  1.1200e-01  2.4000e-02  7.0400e-01  1.2000e-02  1.7380e+00  1.2320e+00
  1.7300e-01  1.0450e+00  6.4100e-01  1.6850e+00  4.4700e-01  2.9870e+00
  3.4350e+00  2.5230e+00  2.1190e+00  4.6420e+00  4.1000e-01  1.6010e+00
  1.5500e+00  2.7800e-01  4.1000e-02  1.1700e-01  3.9960e+00  4.2050e+00
 -1.3400e-01 -5.5300e-01 -1.3400e-01  4.6300e-01  1.1600e-01 -2.5000e-01
  2.5200e-01 -2.7830e+00 -4.3000e-01  3.0000e-02  4.6200e+00  1.0070e+00]
Function at x* is q(x*) = 587.078
Final active set is W* = [1, 5, 7, 9, 11, 3, 13]

n = q_traj.size
x = np.arange(n)
plt.plot(x, q_traj, marker='o',color='black')
plt.xlabel('Step')
plt.xticks(x)
plt.ylabel('q')
plt.title('q optimization per iteration')
plt.grid(True)
plt.show()

Comparamos solución con CVXPY.

# Define and solve the CVXPY problem.
problem_load = loadProblem("lp_blend.mat")
Ae = problem_load["AE"]
n = Ae.shape[1]
Al = -np.identity(n)
A = np.vstack((Ae, Al))
be = problem_load["bE"]
bl = np.zeros(n)
b = np.concatenate((be, bl))
c = problem_load["c"]
G = np.identity(n)
x = cp.Variable(n)
prob = cp.Problem(
    cp.Minimize((1 / 2) * cp.quad_form(x, G) +
                c.T @ x), [Al @ x <= bl, Ae @ x == be]
)
prob.solve()
solution = x.value

print("Method used:",prob.solver_stats.solver_name)
print("Active set optimum value:", round(float(problem_blend.evaluate(xf)),2))
print("CVXPY optimum value:",round(prob.value,2))
print("Normalized euclidean distance:",round(2*np.linalg.norm(xf-x.value)/(np.linalg.norm(xf)+np.linalg.norm(x.value)),2))

[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.
 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
[inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf inf
 inf inf inf inf inf inf]
Method used: OSQP
Active set optimum value: 587.08
CVXPY optimum value: 712.13
Normalized euclidean distance: 0.44

Nótese que CVXPY utiliza el método OSQP. Los valores resultados obtenidos con ambos métodos son similares (distancia normalizada < 0.5). El método del conjunto activo nos da un mejor minimizador.

Observaciones adicionales

Al trabajar con matrices dispersas en optimización con igualdades, pueden surgir problemas de precisión numérica, afectando la estabilidad de la solución, especialmente si las matrices son mal condicionadas o las tolerancias no son adecuadas, lo que puede generar resultados ligeramente diferentes dependiendo del solver y la configuración numérica.

print("x0:",x0)
print("W0:",W0)
print("iteraciones:",x_traj.size)
print("Valor óptimo obtenido:",round(float(problem_blend.evaluate(xf)),2))

x0: [-0.   -0.    0.    0.   -0.   -0.   -0.   -0.    0.    0.   -0.    0.
 -0.   -0.    0.    0.   -0.   -0.    0.   -0.    0.   23.26  5.25 26.32
 21.05 13.45  2.58 10.   10.    0.    0.    0.   -0.   -0.   -0.   -0.
 -0.   -0.   -0.   -0.    0.   -0.   -0.   -0.   -0.    0.    0.    0.
 -0.   -0.   -0.    0.    0.   -0.   -0.   -0.   -0.    0.    0.   -0.
  0.   -0.    0.    0.   -0.   -0.   -0.   -0.    0.    0.   -0.    0.
 -0.   -0.   -0.   -0.    0.    0.    0.    0.   -0.   -0.   -0.   -0.
 -0.   -0.    0.   -0.   -0.   -0.   -0.   -0.   -0.   -0.   -0.   -0.
  0.   -0.   -0.   -0.   -0.    0.    0.    0.    0.    0.   -0.    0.
 -0.    0.    0.   -0.   -0.   -0.  ]
W0: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73]
iteraciones: 228
Valor óptimo obtenido: 587.08

---

VI. Referencias

[1] Nocedal, J., & Wright, S. J. (2006). Numerical Optimization. Springer.
[2] Andreas Wachtel, Curso de Optimización Numérica, Instituto Tecnológico Autónomo de México (ITAM), primavera 2025. [3] Harris, C. R., Millman, K. J., van der Walt, S. J., et al. (2020). Array programming with NumPy. Nature, 585(7825), 357–362. https://doi.org/10.1038/s41586-020-2649-2
[4] Virtanen, P., Gommers, R., Oliphant, T. E., et al. (2020). SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17(3), 261–272. https://doi.org/10.1038/s41592-019-0686-2
[5] Hunter, J. D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering, 9(3), 90–95. https://doi.org/10.1109/MCSE.2007.55
[6] Caswell, T. A., Droettboom, M., Lee, A., et al. (2016). Matplotlib/matplotlib: v1.5.1. Zenodo. https://doi.org/10.5281/zenodo.44579
[7] Diamond, S., & Boyd, S. (2016). CVXPY: A Python-Embedded Modeling Language for Convex Optimization. Journal of Machine Learning Research, 17(83), 1–5. http://jmlr.org/papers/v17/14-541.html
[8] Plotly Technologies Inc. (2015). Collaborative data science. Montreal, QC. https://plotly.com/python/
[9] Python Software Foundation. random — Generate pseudo-random numbers. Documentación oficial de Python 3. https://docs.python.org/3/library/random.html
[10] Andreas Wachtel, Instituto Tecnológico Autónomo de México (ITAM), primavera 2025. Archivo AW_loadProblem.py.